#!/usr/bin/perl
package get_genomic_sequence;

require Exporter;


our @ISA = qw(Exporter);
our $VERSION = 1.00;
our @EXPORT = qw(get_genomic_sequence_coordinates get_genomic_sequence get_genomic_sequence_coordinates_from_index_file);

use strict;
use warnings;
use Carp;
use die_error;
use FileHandle;

sub get_genomic_sequence_coordinates($\%)
{
	use DBI;
	my ($dbh, $genomic_sequence_coordinates) = @_;
	my $array = $dbh->selectall_arrayref(
				"SELECT species, name, seq_beg, seq_len, chrm_offset ".
					"FROM taxon.chromosomes;")
			 or croak DBI::errstr.".";

	for my $line(@$array)
	{
		my ($species, $chrm, $seq_beg, $seq_len, $chm_offset) = @$line;
		$genomic_sequence_coordinates->{$species}{$chrm} = [$seq_beg, $seq_len, $chm_offset];
	}

}





#_________________________________________________________________________________________
#
#       get_genomic_sequence_coordinates_from_panda
# 
# 			get cached indexed genomic data from index file
# 
#_________________________________________________________________________________________
sub get_genomic_sequence_coordinates_from_index_file ($\%)
{
	my ($fh, $genomic_sequence_coordinates) = @_;

	for (<$fh>)
	{
		chomp;
		
		my $line = $_;
		
		#ignore empty or comment lines
		next unless (length($line)  and $line !~ /^#/);

		my ($chrm_description, $acc_start, 
			$acc_len, $seq_start, $seq_len)  = split(/\t/, $line);


		# ignoring ">" in chromosome name
		$chrm_description = substr($chrm_description, 1);
		 
		my @data = split(/:/, $chrm_description);
		my @parts = split(/[: ]/, $chrm_description);
		my $chrm_name = $parts[0];

		#>NT_166479:1:23934:1
		#>Y dna:chromosome chromosome:NCBIM37:Y:1:15902555:1
		my $chrm_offset = 0;
		#if ( @data >= 4)
		#{		
		#	#print $data[-3], $chrm_description;
		#	$chrm_offset = $data[-3] - 1;
		#}		

		$genomic_sequence_coordinates->{$chrm_name} = [$seq_start, $seq_len, 
														$chrm_offset];
	}
}
    


#_________________________________________________________________________________________
#
#       get_genomic_sequence
# 
#_________________________________________________________________________________________
my $cached_fh = undef;
my $cached_filename = "";

sub get_genomic_sequence ($$$\%$)
{
	my ($contig, $beg, $seq_len, $genomic_sequence_coordinates, 
                            $sequence_filename) = @_;
	#print STDERR join(", ", $contig, $beg, $seq_len, $genomic_sequence_coordinates, 
    #                        $sequence_filename), "\n";
    #
    #   use cache
    #


    #
    #   open file and add to cache
    #
    if ($sequence_filename ne $cached_filename)
	{
		close($cached_fh) if defined $cached_fh;
		$cached_filename = $sequence_filename;
		open($cached_fh, "<:mmap", $cached_filename);
	}	

    my $new_contig = undef;
    if (!exists $genomic_sequence_coordinates->{$contig})
	{
		my @parts = split (/\./, $contig);
		my $contig_sans_version = $parts[0];
		#
		# DEBUGG
		if (exists $genomic_sequence_coordinates->{$contig_sans_version})
		{
			$new_contig = $contig = $contig_sans_version;
		}
		else
		{
			#        5|NT_113890.1
			my @contig_sans_chromosome_prefix = split (/\|/, $contig_sans_version);
			if (@contig_sans_chromosome_prefix > 1 &&
				exists $genomic_sequence_coordinates->{$contig_sans_chromosome_prefix[1]})
			{	
				$new_contig = $contig = $contig_sans_chromosome_prefix[1];
			}
			else
			{
				croak ("File $cached_filename has no contig [$contig]");
			}
		}	
	}		

    my ($chrm_beg, $chrm_len, $chrm_offset) = @{$genomic_sequence_coordinates->{$contig}};

    if ($beg < $chrm_offset)
	{
		croak("In contig [$contig], $beg is less than the first " .
										"available sequence position ($chrm_offset).");
	}								
    if ($beg + $seq_len > $chrm_len + $chrm_offset)
	{
		my $end = $chrm_len + $chrm_offset;
		carp ("In contig [$contig], $beg + $seq_len is beyond ".
					"the available sequence length ($end)");
		$seq_len = $chrm_len + $chrm_offset - $beg
	}			

    # 
    # read the file at the specied position
    #
	
    sysseek($cached_fh, $chrm_beg + $beg - $chrm_offset, Fcntl::SEEK_SET);
	my $sequence = "";
	sysread($cached_fh, $sequence, $seq_len);

    if (length($sequence) != $seq_len)
	{
		my $readlen = length($sequence);
		croak("Only able to read ($readlen) of genomic sequence.\n" .
				"Parameters = contig[$contig] [$beg+$seq_len] from [$cached_filename]\n".
				"Contig params : [b = $chrm_beg, l = $chrm_len, offset = $chrm_offset]");
	}

    return $sequence, $new_contig
}

1;
